Authors: J. Coelho, S. Ammer
library(ggplot2)
library(Momocs) # package for outline and curve analysis
## This is Momocs 1.1.0
##
## Attaching package: 'Momocs'
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:base':
##
## table
setwd("/Users/del/ACADEMY/WIP/HumerusProject/Photos_final/") # location of my data
rawLPF <- import_tps(tps.path = "LPF/LPF.TPS", curves = T) # read data Left, Posterior, Females
rawLPM <- import_tps(tps.path = "LPM/LPM.TPS", curves = T) # read data Left, Posterior, Males
For some reason, I can’t use Momocs::Out() directly on my imported TPS. Tried many ways, always obtaining different errors. I thought this might be because the p of coordinates is not consistent. So I resample them. In my understanding coo_sample() works well for open curves, while I preferred coo_interpolate() for my closed outlines, since they had very inconsistent p.
# Rebuild datasets
# For females
for (i in seq_along(rawLPF$cur)){
# open (curves)
rawLPF$cur[[i]][[1]] <- coo_sample(rawLPF$cur[[i]][[1]], n = 16)
# closed (outlines)
rawLPF$cur[[i]][[2]] <- coo_interpolate(rawLPF$cur[[i]][[2]], n = 24)
# keep the IDs
names(rawLPF$cur)[[i]] <- names(rawLPF$cur)[[i]]
}
## For males
for (i in seq_along(rawLPM$cur)){
# open (curves)
rawLPM$cur[[i]][[1]] <- coo_sample(rawLPM$cur[[i]][[1]], n = 16)
# closed (outlines)
rawLPM$cur[[i]][[2]] <- coo_interpolate(rawLPM$cur[[i]][[2]], n = 24)
# keep the IDs
names(rawLPM$cur)[[i]] <- names(rawLPM$cur)[[i]]
}
I tried to work with my curves together but for some reason, I also got errors trying to Out() these. So first I will split my open curve (TE: Throclear Extension) from my closed outlines (OF: Olecranon Fossa), but combine my males and females datasets to see the differences among these groups.
TE.f <- list()
for (i in seq_along(rawLPF$cur)){
TE.f[[i]] <- rawLPF$cur[[i]][[1]]
names(TE.f)[[i]] <- names(rawLPF$cur)[[i]]
}
TE.m <- list()
for (i in seq_along(rawLPM$cur)){
TE.m[[i]] <- rawLPM$cur[[i]][[1]]
names(TE.m)[[i]] <- names(rawLPM$cur)[[i]]
}
Sex <- rep("Female", length(TE.f))
TE.f <- Out(TE.f, fac = as.data.frame(Sex))
Sex <- rep("Male", length(TE.m))
TE.m <- Out(TE.m, fac = as.data.frame(Sex))
TE <- combine(TE.f, TE.m)
###
OF.f <- list()
for (i in seq_along(rawLPF$cur)){
OF.f[[i]] <- rawLPF$cur[[i]][[2]]
names(OF.f)[[i]] <- names(rawLPF$cur)[[i]]
}
OF.m <- list()
for (i in seq_along(rawLPM$cur)){
OF.m[[i]] <- rawLPM$cur[[i]][[2]]
names(OF.m)[[i]] <- names(rawLPM$cur)[[i]]
}
Sex <- rep("Female", length(OF.f))
OF.f <- Out(OF.f, fac = as.data.frame(Sex))
Sex <- rep("Male", length(OF.m))
OF.m <- Out(OF.m, fac = as.data.frame(Sex))
OF <- combine(OF.f, OF.m)
My pics were (anatomically) upside-down. So I use the very useful coo_rotate() by Pi radians (180º). I have done some landmarks-based analysis before were Procrustes was a mandatory step. I am not sure if it is with open curves and closed outlines, however I am also doing a superimposition step here.
TE.aligned <- TE %>% coo_rotate(.,theta = pi) %>% fgProcrustes() %>% Opn()
OF.aligned <- OF %>% coo_rotate(.,theta = pi) %>% fgProcrustes() %>% coo_close()
TE.aligned %>% stack(., title = "Trochlear Extension (H. sapiens)", fac = "Sex")
ftimes <- rep('red', length(grep('Female', TE.aligned$fac$Sex)))
mtimes <- rep('blue', length(grep('Male', TE.aligned$fac$Sex)))
panel(TE.aligned, borders = c(ftimes, mtimes), names = substr(names(TE.aligned), 7, 9))
TE.aligned %>% npoly() %>% PCA() %>% plot(.,"Sex") # I need to read about this
## 'nb.pts' missing and set to: 16
## 'degree' missing and set to: 5
TE.aligned %>% opoly() %>% PCA() %>% plot(.,"Sex") # read about this
## 'nb.pts' missing and set to 16
## 'degree' missing and set to 5
TE.aligned %>% dfourier() %>% PCA() %>% plot(.,"Sex") # strange results; might not be proper method for this dataset.
## 'nb.h' not provided and set to 12
## ===========================================================================
TE.n <- npoly(TE.aligned, nb.pts = 16, degree = 5)
TE.pca <- PCA(TE.n)
TE.pca %>% as_df() %>% ggplot() +
aes(x = PC1, y = PC2, col = Sex) + coord_equal() +
geom_point() + geom_density2d() + theme_light()
# mean shape:
TE.n %>% mshapes() %>% coo_plot()
## no 'fac' provided, returns meanshape
# mean shape per group:
TE.ms <- mshapes(TE.n, 1)
TE_female <- TE.ms$shp$Female %T>% coo_plot(border = "red")
TE_male <- TE.ms$shp$Male %T>% coo_draw(border = "blue")
legend("topright", lwd = 1, col = c("red", "blue"), legend = c("Female", "Male"), cex = 0.8)
title(main = "Throclear Extension - Mean shapes")
# TPS plots
coo_arrows(TE_female, TE_male); title("Deformations")
## Warning in arrows(coo1[s, 1], coo1[s, 2], coo2[s, 1], coo2[s, 2], length =
## length, : zero-length arrow is of indeterminate angle and so skipped
tps_grid(TE_female, TE_male)
tps_arr(TE_female, TE_male)
tps_iso(TE_female, TE_male)
# Manova analysis
MANOVA(TE.pca, 1)
## PC axes 1 to 2 were retained
## Df Hotelling-Lawley approx F num Df den Df Pr(>F)
## fac 1 0.038565 2.8538 2 148 0.0608 .
## Residuals 149
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MANOVA_PW(TE.pca, 1)
## PC axes 1 to 2 were retained
## $stars.tab
## Female Male
## Female .
##
## $summary (see also $manovas)
## Df Pillai approx F num Df den Df Pr(>F)
## Female - Male 1 0.03713 2.854 2 148 0.0608
# LDA / CVA analysis
TE.lda <- LDA(TE.pca, 1)
## 0.99 total variance
## 2 PC retained
TE.lda
## * Leave-one-out cross-validation ($CV.correct): (57.6% - 87/151):
##
## * Class correctness ($CV.ce):
## Female Male
## 0.7125000 0.4225352
##
## * Cross-validation table ($CV.tab):
## classified
## actual Female Male
## Female 57 23
## Male 41 30
plot(TE.lda)
## NULL
plot_CV(TE.lda)
plot_CV2(TE.lda)
OF.aligned %>% stack(.,title = "Olecranon Fossa (H. sapiens)", fac = "Sex")
OF.aligned %>% panel(., fac = "Sex", names = substr(names(OF.aligned), 7, 9))
OF.aligned %>% efourier(norm=TRUE) %>% PCA() %>% plot(.,"Sex")
## 'nb.h' not provided and set to 6 (99% harmonic power)
OF.ef <- efourier(OF.aligned, nb.h = 6)
OF.pca <- PCA(OF.ef)
OF.pca %>% as_df() %>% ggplot() +
aes(x = PC1, y = PC2, col = Sex) + coord_equal() +
geom_point() + geom_density2d() + theme_light()
# mean shape:
OF.ef %>% mshapes() %>% coo_plot()
## no 'fac' provided, returns meanshape
# mean shape per group:
OF.ms <- mshapes(OF.ef, 1)
OF_female <- OF.ms$shp$Female %T>% coo_plot(border = "red")
OF_male <- OF.ms$shp$Male %T>% coo_draw(border = "blue")
legend("topright", lwd = 1, col = c("red", "blue"), legend = c("Female", "Male"), cex = 0.8)
title(main = "Olecranon Fossa - Mean shapes")
# TPS plots
coo_arrows(OF_female, OF_male); title("Deformations")
tps_grid(OF_female, OF_male)
tps_arr(OF_female, OF_male)
tps_iso(OF_female, OF_male, grid = TRUE)
# MANOVA analysis
MANOVA(OF.pca, 1)
## PC axes 1 to 7 were retained
## Df Hotelling-Lawley approx F num Df den Df Pr(>F)
## fac 1 1.8666 38.133 7 143 < 2.2e-16 ***
## Residuals 149
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
MANOVA_PW(OF.pca, 1)
## PC axes 1 to 7 were retained
## $stars.tab
## Female Male
## Female ***
##
## $summary (see also $manovas)
## Df Pillai approx F num Df den Df Pr(>F)
## Female - Male 1 0.6512 38.13 7 143 9.562e-30
# LDA / CVA analysis
OF.lda <- LDA(OF.pca, 1)
## 0.99 total variance
## 7 PC retained
OF.lda
## * Leave-one-out cross-validation ($CV.correct): (92.1% - 139/151):
##
## * Class correctness ($CV.ce):
## Female Male
## 0.9500000 0.8873239
##
## * Cross-validation table ($CV.tab):
## classified
## actual Female Male
## Female 76 4
## Male 8 63
plot(OF.lda)
## NULL
plot_CV(OF.lda)
plot_CV2(OF.lda)